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The position of propagating population fronts fluctuates because of the discreteness of the in- 
dividuals and stochastic character of processes of birth, death and migration. Here we consider a 
Markov model of a population front propagating into a metastable state, and focus on the weak 
noise limit. For typical, small fluctuations the front motion is diffusive, and we calculate the front 
diffusion coefficient. We also determine the probability distribution of rare, large fluctuations of the 
front position and, for a given average front velocity, find the most likely population density profile 
of the front. Implications of the theory for population extinction risk are briefly considered. 
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I. INTRODUCTION 

Population fronts are ubiquitous in nature. They con- 
trol the speed of a multitude of processes ranging from 
epidemic outbursts to invasion of species to biological 
evolution 0, There is a fundamental difference be- 
tween fronts propagating into an unstable state and those 
propagating into a metastable state (a state which is lin- 
early stable but no n- linearly unstable). The former can 
be extremely sensitive to intrinsic noise in the leading 
edge 0-0I- I n particular, for pulled fronts, the front dif- 
fusion coefficient Df, describing wandering of the front 
position around the position, predicted by the determin- 
istic rate equation, is only logarithmically small with N, 
the number of particles in the front region Q. Fronts 
propagating into metastable states are expected to be- 
have "normally" in this respect, and Df ~ 1/N scaling 
was conjectured [H, |f|. Surprisingly, not much is known 
beyond this conjecture. Yet fronts, propagatinginto 
metastable states, appear in many applications 0,H|, 
For example, they describe colonization and extinction 
waves when an Alice effect is present in the form of a 
critical population size needed for establishment of a lo- 
cal population 0, 0] • Motivated by this and other front 
propagation phenomena - in physics, chemistry and biol- 
ogy - we present here a systematic theory of fluctuations 
of the position of a population front propagating into a 
metastable state. 



II. MODEL 

Let a single population of individuals (we will call them 
particles) reside on a one-dimensional lattice of sites la- 
beled by index i. The population size rii at each site 
varies in time as a result of two types of Markov processes 
■ The first type involves stochastic births and deaths 
with rates X(rii) and /z(rij), respectively. The second one 
is random walk of each particle between neighboring sites 
with rate coefficient Dq. 

We assume that the typical population size scales as a 
large parameter, K 3> 1. Then for n, 3> 1 one can write, 



in the leading order [ll|, [l2| : 

\(m) = vK\(qi) and fj,(rii) = vKp(qi), (1) 

where qi = rii/K is the rescaled population size at site i, 
X(li) ~ Pili) ~ 1) an d ^ is a rate coefficient. Assuming 
also a fast diffusion, Dq 3> v, we can use a continuous 
position coordinate x instead of the discrete index i. Ne- 
glecting for a moment the intrinsic noise, we arrive at a 
deterministic reaction-diffusion equation 



dtq = f(q) + d 2 x q , 



(2) 



where f(q) = X(q) — £t(q), D = DqK 2 is the diffusion 
coefficient of the random walk, and h is the lattice con- 
stant. The position coordinate x in Eq. ([2]) is rescaled 
by the characteristic diffusion length Id = (D jv) 1 ! 2 , and 
time is rescaled by 1/V. We will restore the dimensions 
in our main results by recalling that the front velocity is 
measured in the units of vln = (vD) 1 / 2 , and the front 
diffusion coefficient Df is measured in the units of D. 

A (deterministic) population front corresponds to a 
traveling wave solution (TWS) of Eq. ©, q(x,t) = <?(£), 
£ = x — cot, which interpolates between two homoge- 
neous states described by zeros of f(q). Let one of them, 
q = q* > 0, be located at x — > — oo, whereas the other, 
q = 0, at x — > oo. Each of them is linearly stable: 
f'(q*) < and /'(0) < 0. There is also a linearly un- 
stable state q u , so that < q u < g». The deterministic 
TWS go(£) obeys the equation 



q' ' + c Q q' + /(g ) = . 



(3) 



This TWS is unique up to an arbitrary shift in £. The 
formal definition of the metastable state in this work is 
based on the effective potential V(q) = f(q') dq' of the 
deterministic problem 0, H, |1| • The state q = > is 
advancing (c > 0) when V(q*) > 0. In this case q = 
is called the metastable state. For V{q ie ) < the state 
q = q* is retreating (cq < 0), and so it is metastable. 
Finally, for V{q*) = the front is standing, cq = 0, and 
so the two states q = q* and q = coexist. The typical 
front width is on the order of Id- 



2 



Now we go back to the discrete-lattice model and 
describe stochastic dynamics of the population. Intro- 
duce the multivariate probability distribution P(n, t) = 
P(ni,ri2, ■ ■ ■ ,t) = P(h,rii,t), where i = 1,2,... and n 
denotes the vector of all n's not explicitly written, see e.g. 



Ref. [l(|. For the continuous-time Markov processes of 
on-site births and deaths and unbiased random walk be- 
tween neighboring sites, the master equation for P(n, t) 
has the following form: 



d t P(n, t) = { A (™< _ W*- n * - x > l ) + <"( n ' + n * + 1 » *) _ t A ( n ') + Kn t )}P(n, t)} 

+ £>o J! + 1 ) p ( fi ' n '-i + !» »i - !> f ) + + ^C 6 ' n < - + M) ~ 2njP(n, i)} , (4) 



r 



where we have returned to dimensional units. 



III. SMALL FLUCTUATIONS 

Small fluctuations around the deterministic TWS can 
be described by the (approximate) Fokker-Planck equa- 
tion obtained, via a truncated expansion, from Eq. (j4|). 
see e.g. Ref. (131 ): 



(9 



9 



9g. t dq i+ i 



d 



2K V% dq l+ i 



(qi -Qi+i) P 



P 



(5) 



Here <?(g) = A(g) + p,(q) is a rescaled on-site diffusion 
coefficient in the space of population sizes. Instead of 
dealing with the Fokker-Planck equation directly, we will 
consider the Langevin equation to which it is equivalent. 
Reintroducing a continuous description in space, we ob- 
tain H 



d t q(x,t) 



vf(q) + Ddtq+J^0± V (x,t) 



2q(x,t) Dh 



K 



X(x,t) 



(0) 



where rj(x, t) and x{ x i are independent Gaussian noises 
which have zero means and are delta-correlated both in 
x and in t: 

(r](x, t) V (x',t')) = ( X (x, t) X (x', 0) = S(x - x') 5(t - 0- 

(7) 

The first two terms on the right hand side of Eq. ([5]) are 
the same terms as in the deterministic equation ([2]). The 
last two terms describe multiplicative noises coming from 
the on-site birth and death processes (the third term), 
and from the random walk (the fourth term, which has 



the form of a flux). Importantly, the magnitudes of the 
two noise terms are comparable to each other in the front 
region. The flux term does not appear in the Langevin 
equations used in many papers for modeling fluctuations 
of the front position caused by external noise (l5l [l(| . 
In addition, for external noise the analogs of function 
g(q) are chosen ad hoc dEHU, whereas g(q) in Eq. © is 
uniquely determined by the on-site birth and death rates. 

Pechenik and Levine [llj used the Doi-Peliti formalism 
[l8T | to derive the Langevin equation for a microscopic 
model which includes processes A 2A and random 
walk. This model corresponds to propagation into unsta- 
ble state. However, the derivation would hold for prop- 
agation into a metastable state as well. The Langevin 
equation of Ref. [l7| describes the evolution of the sys- 
tem in terms of the dynamics of a sequence of Poissonian 
states A possible advantage of their description is 

that the random walk noise is "absorbed" by the variables 
and does not appear explicitly. An advantage of Eq. 
is that it provides a direct description of the dynamics 
of the system in terms of the (rescaled) population size 
q(x,t). 

Rcscaling the position coordinate x and time t by Id 
and respectively, we can rewrite Eq. ([5]) as 

d t q(x, t) = f(q) + 8%q{x, t) + N~ X / 2 R(x, 1, q) , (8) 
where 



R(x, t, q) = v/fffe) V(x, t) + d x [^2~q X (x, t)} , (9) 

where N = Kl jj /h 3> 1 is the typical number of particles 
in the front region. The smallness of A" 1 / 2 , the effec- 
tive magnitude of the intrinsic noise, enables one to solve 
Eq. ([8]) perturbatively around the deterministic TWS, 
see Refs. [H|,[l6|. One looks for a perturbative solution 



q(tx,t) = q (t;x) + qi(€x,t) , M < q 



(10) 



where £x = x — cot — X(t), and X(t) describes the noise- 
driven displacement of the front position with respect to 
the deterministic TWS. Going over to the variables £x 
and t and linearizing Eq. (jSJ around qo, we obtain 



d t qi = Cqi + Xq' + N~ X ' 2 R{x, t, q ) 



111) 
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where 



91 



co<% + /'(g ) 



(12) 



and we kept the old variable x in the noise term. 

X(t) can be found from the solvability condition of 
Eq. (flTl) . To this end let us consider the eigenmodes 
of the operator C: C(j) n = a n (f) n . The zero mode 
<Ao (£) = QoiOi f° r which a = 0, is the Goldstone mode, 
coming from the translational invariance of the deter- 
ministic TWS. Importantly, the lowest excited mode is 
gapped [HI, . We expand 

3i(6r,*) = $Z &»(*)&»(&:), 

n 

and fix the definition of X{t) by setting bo(t) = 0. Now 
we multiply Eq. (fTTj) by e Co4 </>o(£) and integrate over £ 
from — oo to oo. As the operator e c °^£ is Hermitian, the 
functions 4> n with different n are orthogonal to each other 
with weight function e c °^, so upon normalization 



where 5 nm is the Kronccker's delta. As a result, 



X = - 



^ 1/2 L°l^(9o) 2 e Co « 



(13) 



(14) 



The front diffusion coefficient Df, in the frame moving 
with the velocity Co, is given by the relation 



Df = lim 



/>! f*dh (X(ti)*(t a )> 



2t 



(15) 



Using Eqs. (fl~5|) . (fT4|) and ([9]) we obtain, after some alge- 
bra (see Appendix), 



D f = 



D 



(16) 



where 



so 



J- a L^{5('?o^) 2 3(go) + g [(^ Co «)'] } 



(17) 

That is, in order to calculate the front diffusion coeffi- 
cient one only needs two ingredients: the deterministic 
TWS <Zo(£) and the function o(g) = X(q) + /2(g) evalu- 
ated on <?o(0- As expected d,|g|, Df /D scales as N^ 1 , in 
stark contrast with the propagation into an unstable state 
where Df/D ~ (lnA^) -3 [7|. The first and second terms 
in the denominator of Eq. (|17[) come from the birth-death 
noise and random walk noise, respectively. We note that 
the numerator and the first term in the denominator ap- 
peared in the calculations of the front diffusion coefficient 
for the case of external fluctuations [HI, . In that case, 
however, g(q) is unrelated to the on-site birth and death 
rates. Importantly, for /'(<?*) < and /'(0) < 0, all the 
integrals entering Eq. (|17[) converge. 



IV. WENTZEL-KRAMERS-BRILLOUIN (WKB) 
APPROXIMATION 



Large fluctuations of the front position render the 
Fokkcr-Planck approximation ([5]) and the ensuing 
Langcvin equation © inapplicable. However, one can 
still exploit the small parameter 1/N via a WKB ap- 
proximation, see e.g. (l3l. I20I [2lj . by making the ansatz 
P(n, t) = exp[— KS(q, t)] in the master equation In 
the leading order in 1/ K, this brings about a Hamilton- 
Jacobi equation dtS + H(q, 3 q 5) = 0. In the continuum 
limit in x, the Hamiltonian functional H = v J dxw. 
Here 



= H a {q,p) - d x q d x p + q (d x p) 



(18) 



is the rescaled Hamiltonian density, x is rescaled by lo, 
and 



H (q,p) = X(q) {e p -l) + Jx{q) (e^ - l) 



(19) 



is the rescaled on-site Hamiltonian [13|,l20(. The last term 
in Eq. (|18[) comes from the random walk noise. Note the 
identities H o {q,0) = 0, d p H (q,0) = f(q), d qp H (q,Q) = 
f(q) and d pp H o (q,0) = g(q). 

Instead of directly solving the Hamilton- Jacobi equa- 
tion for S"(q, t), consider the Hamilton's equations of mo- 
tion 



d t q = d p H Q (q,p) + dlq - 2d x (qd x p) 



d t p 



-d q H (q,p) - dip - {d x pY 



(20) 
(21) 



for the population size g(x, t) and effective momentum 
p(x,t) = hSS/Sq. Deterministic dynamics take place in 
the invariant hypcrplanc p = 0, where Eq. (|20p reduces to 
Eq. ©. Fluctuations involve a non-zero p(x, t), and this 
is the case of our interest. To have a well-defined problem 
for Eqs. (f2"0")l and (|2ip one needs to specify the boundary 
conditions, both in space and in time. The boundary con- 
ditions in space correspond to zero-energy fixed points of 
the on-site Hamiltonian: {q(— oo, t),p(— oo, t)} = {g*,0} 
and {q(oo,t),p(oo,t)} — {0,0}. The boundary condi- 
tions in time involve, in general, specifying kink-like pop- 
ulation size profiles qi(x) and q2(x) at t = and t = r, 
respectively. After having solved the Hamilton's equa- 
tions, one can calculate the action along the trajectory 
of the system in the functional phase space q,p: 



KS = N 



dx 



dt \p[x, t)dtq — w] 



(22) 



We are interested in the probability of observing, in the 
frame moving with the velocity Co, a front displacement 
X during a time interval which is much longer than the 
front relaxation time: ft » 1. Our crucial assumption 
is that this probability is described, in the leading order 
of \/(vt) <C 1, by a traveling wave solution (TWS) of 
Eqs. (|2U|) and (|2"Tj) . This leads to two coupled ordinary 
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differential equations for q(£) and p(£), where £ = x — ct 
and c = X/t is the average front velocity during time r: 



P 



2( qP y 
+ (P) 2 



cq' + d p H (q,p) = 0, (23) 
cp' + d q H a (q,p) = 0. (24) 



The boundary conditions are q(—oo) = g» and q(oo) = 
p(-oo) = p(oo) = 0. As d t q = -af(£) for a TWS, 
Eq. ((221) reduces to 



(25) 



Remarkably, Eqs. ([23]) and ([24]) possess a conservation 
law: 



#«S(c) = -Nut / d£ [cp(0?'(0 + to] 



H a [q(0,p(0] + V'P' - q(p') 2 =C = const 



(26) 



where, by virtue of the boundary conditions, C = 0. 

Although especially suitable for dealing with large fluc- 
tuations of the front position, the WKB theory also cap- 
tures typical, small fluctuations. To see it, let us linearize 
Eqs. (J2S1) and ([24]) around the deterministic TWS ©: 

9(O=90(O+«(O and c = c + <5c, (27) 
where |u| ~ \p\ ~ \Sc\ <C 1. We obtain 

Cu = -Scq' -g(q Q )p + 2(q Q p'y , (28) 
£*p = (29) 

with zero boundary conditions at ±oo. The operator 
£ was defined in Eq. ([T2"]). whereas the operator £* = 
d| — coc% + /'(go) is adjoint to £. 

For the small deviations, the action ([2~5]) can be ex- 
panded to second order. The linear terms cancel out, 
whereas the quadratic terms yield, after some algebra, 



/oo 
-OO 



2#(<7o)p 2 



qo(p' 



(30) 



As one can see, the knowledge of p(£) alone is sufficient 
for evaluating the action (and the front diffusion coeffi- 
cient, see below) for small fluctuations. Therefore, let us 
consider Eq. ([2T)]) . One independent solution of this equa- 
tion, pi(£) = e c °^g (£), obeys the zero boundary condi- 
tions at ±oo [22 ]. Therefore, p(£) = poPi(0 with an a 
priori unknown constant pq. To determine this constant, 
we turn to Eq. ([28]) for u(£), where p(£) enters the effec- 
tive forcing term 

F(0 = -6cq' Q (0 -Po5[9o(0]Pi(0 + 2po[go(OPi(0]'- 

Demanding zero boundary conditions at ±oo, one can 
obtain a solvability condition: /_ <i£ e co ^<? (£)-F(£) = 0. 
This yields 



Po 
5c 



II 



1°^ d!i{(q' e^)2g(q ) + 2q [(e^)']' } 



Then, fromEq. KS = (1/4)s Nvt(5c) 2 , where s is 
defined in Eq. (]17p . and 5c is still dimensionless. Recall- 
ing that 5c = X/t — cq, we see that the probability distri- 
bution P(X, t) is a Gaussian in the frame moving with 
the deterministic front velocity Co- Normalizing this dis- 
tribution to unity we obtain, already in the dimensional 
variables, 



P(X,r)~ (^Ly^exp, 



s NX< 



(32) 



where \X/t - c \ < (vD) 1 / 2 . Equation ((35J) coincides 
with the Green function of a diffusion equation, so the 
front diffusion coefficient Df = D/(sqN), in a perfect 
agreement with Eq. (|16[) . derived from the Langevin 
equation. 

Rare large fluctuations of the front position correspond 
to \p\ > 1 in the WKB formalism. This demands a full 
nonlinear solution of Eqs. ([2"3"]l and ([M]) which in general 
is only accessible numerically. Examples of such calcula- 
tions are presented in the next section. 

We should emphasize that our results on fluctuating 
extinction fronts are conditional on the absence of at- 
tempted extinction events in the region to the left of 
the front. This conditioning is especially important for 
Co < 0. Here an extinction front appears already in de- 
terministic description, and formation of critical nuclei 
of the state q = is possible || [l3| . 



V. EXAMPLE 

Consider three on-site reactions: ^4 — !► 0, 2A — > 3A 
and 3 A — > 2 A, with rate coefficients p,o, Xq and ao, 
respectively [12, H3, HE H3]-_ It is convenient to per- 
form the rescaling so that X(q) = 2q 2 and p(q) = 
7<7 + q 3 . where v = 3Aq/(8oo), K = 3Ao/(2oo), and 
7 = 8/iocro/(3A§). The on-site dynamics exhibits bista- 
bility at 5 2 = 1 — 7 > 0. In this case the zeros of function 
f{q) = —q {q — q u )(q — q*) describe two stable fixed points 
of deterministic theory, and = 1 + 5, and an unsta- 
ble fixed point q u = 1 — 5 such that < q u < q*. The 
deterministic TWS q (x - c t) is [J, HQ 



Qo (0 = 



1 



1 + e (l+5)«/V2 



CO 



3(5-1 
~V2~ 



(33) 



(31) 



so the deterministic front is advancing at 1/3 < 5 < 1, 
retreating at < 5 < 1/3 and stationary at 5 = 1/3. The 
function g(q) = X(q) + p(q) = j q + 2q 2 + q 3 . 



A. Small fluctuations 

Within theory of small fluctuations, one only needs to 
calculate sq = sq(5) from Eq. ([T7]) which enters the front 
diffusion coefficient Df — D/(sqN). All the integrals 
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can be evaluated analytically. The result is a bit cum- 
bersome, however, so we only present it graphically, see 
Fig. [U and give its asymptotes: 



so (5) 



2V25 
9 

5yg _ 
6 



178\/2a 2 
135 



95V2 (1-8) 
72 



0< <5« 1, 

o<i-<y«i. 



(34) 



The latter asymptote corresponds to (j,q<jo -C Aq. For 
(5 = 1/3 (when the deterministic front is stationary) we 
obtain so = a/2/6. 

As one can see from Fig. [TJ the front diffusion Df = 
D/(sN) for the extinction fronts is stronger than for 
the colonization fronts. Furthermore, Df diverges as 
1/6 when 6 — > 0. This happens because, in this case, 
the integral in the denominator of Eq. (|17p diverges at 
£ = — 00. The divergence comes from the first term of 
the integrand, describing the birth-death noise. Notice 
that 6 = is the (saddle-node) bifurcation point of this 
model, and it is hardly surprising that fluctuations are 
enhanced in the vicinity of this bifurcation. Actually, 
very close to the bifurcation point our perturbation the- 
ory breaks down. Indeed, it is only valid when the typical 
value of X is much less than the deterministic front ve- 
locity ~ (vD) 1 / 2 . By virtue of Eqs. (0 and ©, this 
boils down to a strong inequality N5 S> 1. 

On the other hand, Df behaves regularly at 6 = 1, 
that is when the reaction A — s- is absent. Here Sq takes 
its maximum value 5-\/2/6, so Df reaches its minimum 
for this model. Notice that, for the reactions 2A 3A 
the state q = is marginally stable, but unstable with 
respect to any finite perturbation q > 0. 
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FIG. 1: (Color online) Function so, defined in Eq. (|17[) . versus 
S for processes A — >• 0, 2A <=^ 3A, and random walk. The solid 
and dashed segments distinguish between the colonization and 
extinction fronts, respectively. 



B. Large fluctuations 

Now consider large deviations of the front position. We 
solved nonlinear equation s (l23l) and (|24|) for the TWS 
by the shooting method [25|. For given 5 and c we 
first establish analytically, by linearization, the asymp- 
totic behavior of the solutions at large negative £, where 



I q — q*\ <C 1 and \p\ 1. This leaves us with two arbi- 
trary constants instead of four. By fixing one of them, 
we fix the location of the solution in £, and are therefore 
left with a single shooting parameter. The conservation 
law ([26)) is used for accuracy control. 
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FIG. 2: (Color online) Solid lines: the most probable pop- 
ulation size profiles for rescaled front velocities c = 2v2/3 
(a) and c = — y2/3 (b). Parameter 5 = 1/3, so that the 
deterministic front (dashed line) does not move. 

Figure [5] which shows the numerically found pro- 
file for 5 = 1/3 and c = 2^2/3 (a) and c = -y/2/3 (b). 
Recall that the deterministic front is stationary in this 
case, so such high front velocities demand large fluctua- 
tions. The profile - the most probable population 
size profile for this front velocity - is steeper than the 
deterministic profile in Fig. [5^, and less steep that the 
deterministic profile in Fig. [2}d. 

Figure [3] shows the numerically found rescaled "action 
accumulation rate" 



ds 
di 



d£(cpq' + w) 



(35) 



[see Eq. (|23|) ] versus c for 5 = 1/3. For small |c| the 
numerical results agree with the linear theory prediction 
ds/dt = V2(Sc) 2 /2A. The tails of the distribution ~ 
exp[— KS(c)] = cxp(— Nurds /dt), however, are strongly 
non-Gaussian: under-populated for fluctuations pushing 
the front to the right (toward colonization), and over- 
populated for fluctuations pushing the front to the left 
(toward population extinction). That is, there is a much 
higher probability, at cq = 0, to observe an unusually 
fast extinction front, than an unusually fast colonization 
front. This feature, however, can be model-dependent. 
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FIG. 3: (Color online) Symbols: numerically found "action 
accumulation rate" ds/dt, see Eq. (|35[) . as a function of the 
rescaled front velocity for S = 1/3. Dashed line: the predic- 
tion of linear theory, ds/dt = \/2 (<5c) 2 /24. 



mediately obtain pi = 1 — Xo/i and p r = xq/L. For 
the extinction via the left boundary, 7} coincides, in the 
leading order, with the average first-passage time. The 
latter. 




is well known, see e.g. Ref. [27|, P- 47. Typically, 
T; ~ L 2 /Df in this case. For the extinction via the right 
boundary, the conditional mean time to extinction T r is 
expected to be exponentially long: much longer than the 
first-passage time for that boundary. 



VII. SUMMARY 



VI. IMPLICATIONS FOR POPULATION 
EXTINCTION 

Now let the spatial system size be large but finite: 
L Id, with impenetrable (reflecting) boundaries at 
x = and x = L. We still assume bistability, with 
q = q* and q = being two linearly stable fixed points 
of the deterministic description. The zero population 
size throughout the system is an absorbing state of the 
stochastic dynamics and, in the absence of immigration, 
the population ultimately reaches this state with proba- 
bility one. 

Let at t = a population front be located at x = xq, 
where < xq < L and not too close to the boundaries. 
As before, we assume that the populated state q = q* is 
on the left. The mean time to extinction, T can be very 
different depending on whether the deterministic front is 
retreating, cq < 0, or advancing, Co > 0. For cq < 0, the 
population goes extinct essentially deterministically, and 
T ~ aco / 1 Co | scales as the system size L. For c > the 
front first populates the whole system, around q = q*, 
with probability close to 1. This happens on a time scale 
(L — xo)/co. Then the probability that the population 
goes extinct by time t grows with t very slowly, as T is 
expected to be exponentially long (ill . [2G | . 

Now, what happens in the special case of c = 26] ? 
Here the front position wanders because of the intrin- 
sic noise. As a result, there are two routes to extinc- 
tion. In the first route the wandering front first reaches 
the left boundary x — where the population gets ex- 
tinct rapidly. In the second route, the front first reaches 
x = L. Here the population develops a long-lived quasi- 
stationary distribution around q = q*, and the probabil- 
ity togo extinct by time t grows with time very slowly 

mm 

The probabilities pi and p r of each of the two routes 
to occur, and the corresponding conditional mean times 
to extinction 7] and T r are described, in the leading or- 
der, by a classical first passa ge p roblem for random walk 
(of the front) in an interval jlOl. l27j|. Therefore, we im- 



We have developed a systematic theory of velocity 
fluctuations of a population front propagating into a 
metastable state. These fluctuations result from the in- 
trinsic noises of reactions (births, deaths, etc.) and ran- 
dom walk, and we have identified the contributions of 
each of these noises to the velocity fluctuations. For typ- 
ical, small fluctuations we have derived, from the un- 
derlying microscopic model, the Fokker-Planck equation 
and the effective Langevin equation. We have calculated 
the front diffusion coefficient using two different meth- 
ods. The first method involved a pcrturbativc solution of 
the Langevin equation in the weak-noise-limit, 1/N <C 1. 
The second method employed a perturbative solution of 
WKB equations in the form of a traveling front which 
involves both the rescaled population size q(x, t) and the 
effective momentum p(x, t). As expected, the front diffu- 
sion coefficient scales with N as 1/N. We have also eval- 
uated, with a WKB theory, the probability of rare large 
fluctuations of the front. Here the probability distribu- 
tion of the front position, in the frame moving with the 
deterministic front velocity, turns out to be non-Gaussian 
and strongly asymmetric. We have also considered some 
implications of the theory for evaluating the population 
extinction risk. Our approach is readily extendable to 
multiple-step on-site processes such as nA — > (n + k)A , 
where n and k are integers, and k is not too large [20ll23j|. 
Finally, most of our results (one notable exception be- 
ing the extinction aspects considered in Sec. VI) remain 
valid when both of the two linearly stable fixed points of 
f(q) describe populated states, as in the spruce budworm 
model, see e.g. Ref. p. 7. 
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Appendix 

Here we give some detail on the algebra leading to 
Eq. ([To]). Substituting Eq. (JUJl in Eq. ([T5|). wc obtain 



Df = lim 

J t-KX) 



Ho Ho d & Jo dh J* dt 2 g (£i)9o(6)e Co( « 1+ ^ (R{x 1 ,h,q )R(x 2 ,t 2 , go)) 



I 



(Al) 



Now we perform the ensemble averaging using Eq. 
As 77 and \ are uncorrelated, we are left with two terms: 



Df = Dfi + Df2- The first term comes from the intrinsic 
noise of births and deaths: 



Dfi = lim 



J-oo <%i J-00 «2 Jo dh Jo dt 2 q'oiZM&e^^ v/SWI% - - 



2iVt 



I 



(A2) 



which yields 



Dfi = 



JZoWoe^gfo) 



2N 



Now we calculate the second term, coming from the ran- 
dom walk noise: 



(A3) 



D 



/2 



lim 



He d^i Ho ^2 Jo d h Jo d h 9 (a)9o(6)e Co( « 1+ « 2) <^ 1 V^O^TxCfi,*!) % y^(6) X (£ 2 ,t 2 ) ) 



2At 



We have 

(% fv^(6)x(6,«i)l % fV2®(6)x(6,*2)l ) 



V?o(a)<?o(6) (x(a, *i)x(6, *2)> 



= 2*(f 1 -t 2 )^ li& [90(6) 5(6 -6)] 

= -2<5(t! - ta) % [90(6) -5' (6 - 6)] • (A5) 

The double integration over time of 5(£i — £2) yields t 
which cancels with t in the denominator. Integration 
over £2 gives 



^2<Zo(6)e CoS2 <S'(a-6) 



-2K(aK° 51 ]'- 

(A6) 



Now we calculate the integral over £i : 



dCi^iK * 1 {«>(&) [?o(6)e c ^] ' }' . (A7) 



Integrating by parts, wc can rewrite this integral as 
2JZ d^qo(0 He cai )f- Therefore, 



D 



J^qoUqoe^)' 



/2 



iV 



(A8) 



Summing up the contributions (| A3|) and (|A8|) . we obtain 
D/ = D/(sqN), where so is defined in Eq. (fT7|). and we 
have restored dimensions by reintroducing D. 
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